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Abstract 

A theoretical model for stratified epithelium is presented. The viscoelastic properties of the tissue 
is assumed to be dependent on the spatial distribution of proliferative and differentiated cells. 
Based on this assumption, a hydrodynamic description for tissue dynamics at long-wavelength, 
long-time limit is developed, and the analysis reveals important insight for the dynamics of an 
epithelium close to its steady state. When the proliferative cells occupy a thin region close to the 
basal membrane, the relaxation rate towards the steady state is enhanced by cell division and cell 
apoptosis. On the other hand, when the region where proliferative cells reside becomes sufficiently 
thick, a flow induced by cell apoptosis close to the apical surface could enhance small perturbations. 
This destabilizing mechanism is general for continuous self-renewal multi-layered tissues, it could 
be related to the origin of certain tissue morphology and developing pattern. 
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I. INTRODUCTION 


Biological tissues are viscoelastic. Concepts devel^ed for describing passive viscoelastic 
materials have been applied to biological tissues [l, 2|- However, active processes such as 
forces generated by the cells, cell movement, and cell proliferation make biological tissues 
different from passive viscoelastic materials. Thus biological tissues fall into the category of 


active soft matters [3[. On the other hand, junctions between neighboring cells prevent free 
diffusion of cells in a tissue. It makes biological tissues different from other active biological 
systems such as bacterial colonies and cell cytoplasm. As a result, biological tissues become 
interesting systems to study complex collective motion and morphogenesis. For example, 
it is known that the homeostasis state of a tissue often has undulating surface to increase 
the contact area with its environment. Understanding the formation and maintenance of 
such spatial structure is thus important for fundamental research of active matters and for 
applied research such as fabricating functional surfaces and flexible electronics j4|. 

Experiments have demonstrated that at short time scales a tissue is solid-like due to the 
adhesion junctions between cells jsj. Since the residual stresses can relax due to turnover of 


B. 


junction proteins [ 2 |, 1 ^ and cell rearrangement induced by cell division and cell apoptosis 
a tissue is liquid-like at long time scales [l, 2, 8 -ll|. This suggests a characteristic time r 
such that at time scales smaller than r a tissue behaves like a solid with shear modulus E. 
On the other hand, at time scales larger than r, a tissue behaves like a fluid with viscosity 

^ ^ BQ, 


77 . Typically E ~ 10^ — 10^ Pa [lOl], and r ~ 10 — 10^ s 
order of 10^ — 10® Pa s [I 


therefore rj = Et is on the 


In this study we focus on stratihed epithelium, a multi-layered continuous self-renewal 
tissue 13|. Epithelium often forms the outermost layer of skin and mucous membrane, and 
acts as a barrier separating the outside and inside of a multicellular organism |l^. On the 
top of a stratihed epithelium is a free apical surface, on the bottom is a basement membrane 
attached to the underlying connective tissue composed primarily of collagen and elastin 
hbers. Typically the thickness of stratihed epithelium is ~ 100 pm and the thickness of 
the underlying connective tissue is ~ 1 mm 1^. The proliferative cells are often localized 


near the basa’ 
cell niche 


lamina, suggesting the existence of a special microenvironment called stem 
ITj . Above the proliferative cells are terminally diherentiated cells. They are 


functional cells that do not divide, instead they undergo programed cell death (apoptosis). 
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Intuitively, different cells should have different mechanical properties. As supported by 
recent experiments, stem cells and tumor cells often have smaller stiffness than normal 
differentiated cells [l8|, ll9|. This suggests that the mechanical properties of a tissue should 
also depend on local tissue composition j?]. The effects of such inhomogeneity on tissue 


dynamics is the focus of this study. 

Undulating/hngering structures that are important for biological functions are often seen 
in the apical surface or basal membrane of a stratihed epithelium. For example, skin wrinkles 
or folds are formed when a skin is deformed by external/residual stresses with wavelength 


1 mm 


3 , 20 


22| . At length scale ~ 100 pm, rete pegs interdigitate with dermal papillae 


are observed in the interface between skin epithelium and stroma 
have also been observed in the epithelium of mucous membrane 


23 


20 


25l . Similar structures 


26 


Studies have 


shown that elastic forces could be important during the formation of the aformentioned 
structures. By modeling a stratihed epithelium as an elastic solid in contact with the un¬ 
derlying elastic connective tissue, the formation of skin wrinkles was associated with the 
buckling instability originated from the competition between bending energy of epithelium 


and stretching energy of connective tissue 


15 


2l|. The rete peg structure was suggested to 


arise from the incompatible growth of epithelium and connective tissue 29|. Similar mech 


anism can also explain the formation of hngerprint 


30| . crypt |3l|, and surface wrinkling 


of tumor spheroids 32|. On the other hand, hydrodynamic instabilities are also likely to 
produce undulating morphologies in stratihed epithelium. For example the formation of rete 


peg structure could arise from proliferation-induced stresses 


33| . and an essential condition 


of this result is that the net cell proliferative rate deceases with the distance from basal 


lamina. Some models considered more details such as cell lineage dynamics 


34| . viscoelas¬ 


ticity of the underlying stroma 35|], and the ehect of growth factors inside a tissue 3J, l35 |. 
Overall, several interesting structures and patterns can be understood by a continuum de¬ 
scription of tissues without describing the movement of individual cells and biomolecular 
events in details. The mechanical stresses induced by cell turnover inside a tissue 7|, ll2|] is 
sufficient to do this job. 

The aim of this work is to study theoretically how inhomogeneity of mechanical properties 
in a tissue ahects the tissue dynamics. For illustrative purpose we consider a stratihed 
epithelium on a rigid stroma. We begin with a simple model in which the mechanical 
properties of the tissue is assumed to be homogeneous. With proliferative cells occupying 
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the lower part, and terminal differentiated cells occupying the upper part of the tissue, this 
model has a homeostasis state with a flat apical surface. The linearized dynamics close to 
the homeostasis state in this model already reveals interesting effects from flows and cell 
proliferation. For example, cell division/apoptosis speed up tissue relaxation towards the 
homeostasis state; but this effect is relatively weak when the wavelength of the perturbation 
is of the order of the thickness of the tissue. 


Due to frequent cell division events, tissue viscosity in the region rich in proliferative 
cells should be different from the rest of the tissue. This feature is included in our second 
model. Such a modification indeed leads to a different dispersion relation for tissue dynamics 
close to the homeostasis state. Especially, when the proliferative region becomes sufficiently 
thick, cell apoptosis close to the apical region could induce a flow that strongly hinders the 
relaxation towards the steady state. This effect is especially significant for perturbations with 
wavelength comparable to tissue thickness (~ 100 pm for typical stratihed epithlium). Since 
our model is quite general, we expect this destabilizing mechanism to exist in continuous 
self-renewal tissue, but it is not easy to be observed in normal healthy tissue because the 
proliferative cells only occupy a thin region. However, during development or regeneration 


there is often a transient increase of proliferative cell inside tissues 


36l |. As a result, the new 


effect described by our second model could happen, and it is even possible that undulating 
apical surface could be generated by this mechanism. 


This article is organized as follows. Both models are presented in Sec. [TTl where we 
also discuss the homeostasis states and linearized dynamics close to the steady states. In 
Sec, mil we discuss the flow in an epithelium that is induced by cell division/apoptosis. Using 
the result derived in the Appendix, we also discuss the relation between the dimensionless 
parameter r = ro/rs and rjrei = Vd/vs, where and rs are cell apoptosis rate and cell 
division rate, and rj]:,, rjs are viscosity in region rich in differentiated cells and proliferative 
cells, respectively. In Sec. HVl based on these results, we make a conclusion and propose 
possible experimental tests to validate our theory. The Appendix presents the relaxation 
rate in our second model, then a derivation of the tissue viscosity from a hydrodynamic 
model is presented. 
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II. MODELS 


A. First model: stratified epithelium with constant viscosity 


Consider a basel lamina sitting on the a:|/-plane. An epithelinm grows from this snrface 
into z > 0 region. On top of the epithelium is a free apical surface. We are interested 
in the steady state of this system and the linearized dynamics close to this steady state. 


Since epithelial growth typically takes several days to complete 


23 


24j |. in this article we 


focus on the long time behavior of our model. Therefore processes with relaxation time not 


longer than the time scale r introduced in Sec. [I] have all decayed aw^. 


’he elasticity of 


the tissue can be neglected, and the tissue behaves as a viscous fluid j8l-lll|. The viscosity 
of the tissue depends on cell-turnover events and the rearrangement of junction proteins 
{ 2 , 6|. Experimentally, the viscosity of a tissue can be estimated from the elastic constant 
E and relaxation time r through the relation rj = Et. 

Since the dynamics is slow, the inertia of the tissue is also negligible. Imposing incompre- 
sibility for simplicity, the equation of continuity, force balance, and consitutitive equation 
are 

divi = kp, ( 1 ) 


diO'ik 0 , 


( 2 ) 


^ik P^ik T ‘^PVik- 


(3) 


Here Vi is the i-th component of the flow velocity field, Vik = {diVk + dkVi)/2 is the strain 
rate tensor in the tissue, p is the tissue pressure, p is the viscosity of the tissue, and kp is the 
(net) cell proliferation rate in the tissue. In general kp is a decreasing function of 2 because 
both nutrient and proliferative cells are relatively abundant near the basal lamina. In this 
model p is chosen to be a constant for simplicity. However, in general p should depend on 
local tissue composition through forces generated during cell division and cell apoptosis j^, 
this will be considered in Sec. Ill HI when we introduce our second model. 

We focus on the situation when the tissue is soft compared to the stroma, therefore the 
stroma is modeled as a rigid substrate for simplicity, and the flow field vanishes at z = 0. 
More general model that includes the viscoelastic properties of the stroma will be deferred to 
future studies. The apical surface of the tissue is assumed to be free, and the free boundary 
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condition that includes contribution from surface tension of the tissue and external pressure 
is applied. 

The steady state velocity held of Eqs. 0-0 with these boundary conditions is u* = 
Vy = 0 and v*{z) = kp{z')dz'. The superscript * is used to denote the steady state values. 
Since v* should vanish at the apical surface, we have the following steady state condition 

pH* 

0 = / kp{z)dz, (4) 

Jo 

where H* is the thickness of the tissue in the steady state. Since kp is a decreasing function 
of z, for the solution of Eq. (jl]) to exist, the tissue needs to have kp{0) > 0 and kp{H*) < 0. 
Intuitively, this means that cell division events in the lower part of the tissue should be 
balanced by cell apoptosis in the upper part of the tissue. We emphasize that our model 
is independent of the details of kp{z), the only constraint for kp{z) is to have steady-state 
tissue height H* > 0 which satishes Eq. (j3]). 

A homeostasis state of a tissue is a steady state that is stable against small perturbations. 
To check if the steady state in this model corresponds to a homeostasis state, we consider 
a perturbation that takes the tissue height from H* to H{x,t) = H* -|- 6H{x,t). Similarly, 
Vk = 'VkJ~ p = p* + dp, etc. The linearized equations for the perturbed system are 

didvi = 0, (5) 

and 

- dkSp + rjdididvk = 0 . ( 6 ) 

Note that the perturbed tissue is translationally invariant in the ^/-direction. 

For a perturbation with dH, dvi, dp ~ e*<? 3 ;+a;t no-slip boundary condition at z = 0 

gives dvx\z=o = dvz\z=o = 0. The linearized boundary conditions at z = H* are 

- dp\^=H* + 2ridMz\z=H* = adxdxdH, (7) 

and 

V {dJvx\z=H* + dxdvz\z=H*) = -2r]kp{H*)dxdH, ( 8 ) 

where a is the surface tension of the apical surface. The above two equations describe the 
stress balance in the normal and tangential directions of the apical surface. There is also a 
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linearized kinematic boundary condition connecting the evolution of the tissue surface and 
the flow field 

dt6H = 5 v,1=h* + K{H*)6H. (9) 

The dispersion relation u:{q) can be calculated from the linearized equations and boundary 
conditions of the perturbed system, straightforward algebra gives 


^(*?) ^mechiji) '^^phyiSl)- 


( 10 ) 


The first term 


2qH* — sinh 2qH* 


^ 2{qH*y + (1 + cosh 2qH*) ^ ^ 

is the same as the dispersion relation for the surface of a low Reynolds number simple fluid. 
In the limit of large qH* this term approaches aq/2ri |37|. As Fig. [1] (a) shows, Umech is 
negative for all q, this means that the surface tension of the apical surface tends to stabilize 
the steady state. The second term is 

1 + cosh.{2qH*) 


^phy — 




( 12 ) 


2{qH*y + 1 + cosh{2qH*) 

It is always negative since kp[H*) < 0, i.e., the apical region of the tissue is rich in terminal 
differentiated cells which do not divide but undergo apoptosis. The meaning of this term is 
intuitive: cell proliferation and cell apoptosis help to bring a perturbed tissue back to the 
steady state. Note that for given kp{H*), ujphy has a minimum at wavenumber q = E/H*, 
where S is the solution of coshS = SsinhS. Numerically one hnds that S = 1.2 and 
Uphy{q = E/H*) = 0.7kp{H*). 

In this simple model, it is convenient to think that the relaxation of the tissue towards the 
steady state is due to apical surface tension and cell proliferation. Contribution from surface 
tension grows as the wavenumber q of the perturbation increases, approaching aq/2ri in the 
limit of large qH*. Contribution from cell proliferation has a minimum at qH* 0{1). 


B. Second model: effect of viscosity heterogeneity 


In a typical stratified epithelium, proliferative cells are located close to the basal re- 
3, 12]. Due to their small abundancy, typically the division rate of proliferative cells 

Since stress relaxation 

Q. 


gion 

is large compared to the apoptosis rate of differentiated cells 


in a tissue strongly depends on the division and apoptosis of the cells [7[, the viscosity in 
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regions rich in proliferative cells should be smaller than the viscosity in the other regions 
of the tissue. Detailed discussion of the relation between cell division/apoptosis and tissue 
viscosity is presented in Sec. IIIIBl and the Appendix. 

To consider position-dependent viscosity we modify our model in Sec. Ill Al bv assuming 
that all the proliferative cells reside in the region 2 ; < hs, and differentiated cells are located 
at 2 ; > hs- This is clearly an over simplihcation, nevertheless it captures the basic tissue 
composition in a stratihed epithelium. The net cell proliferation rate inside the tissue is 
simply kp = rs for z < hs and kp = —r^ for z > hs- Here rs is the cell division rate and ro 
is the cell apoptosis rate. From Eq. (|1]), the proportion of proliferative region in the steady 
state is related to rs and by 


hs ^ f 
H* 1 + r' 


(13) 


where r = r^/rs- We emphasize that the control parameter in typical experiments is r . 
For example, r can be increased by increasing the apoptosis rate of differentiated cells 


For given r and hs, we can use Eq. flT^ to obtain the steady state tissue height H*- As a 
reminder, note that in this model kp{H*) = — 

With the aforementioned simplihcation, in this model we choose f] = rjs for z < hs and 
r] = rjD for 2 ; > hs- The geometry and boundary conditions for the linearized tissue evolution 
equations are the same as those for our hrst model, except for additional velocity and stress 
continuity conditions at z = hs- The steady state of our second model therefore becomes 
V* = V* = 0, and n* = rsz for < hs, v* = rshs — r^z for hs < z < H*- The continuity 
of normal stress at hs leads to a discontinuity of tissue pressure in our model that should 
be smoothed out in a real tissue because the boundary between proliferative region and 
differentiated region is diffuse, not inhnitely sharp. 

On the other hand, in our second model the dynamics close to the steady state are 
qualitatively different from that of the hrst model. Straightforward calculation for a small 
perturbation of the form 5H ~ the steady state of our second model gives 

the relaxation rate cu. We hnd that u can still be expressed as the sum of two terms oJmech 
and bJphy, but these two terms are diherent from the counterparts of our hrst model. 


^mech Kmechi^hH ,r,T]j-gij 


2vd 


(14) 


^phy I^phyiSlH 


( 16 ) 









The explicit expression of the dimensionless functions Kmech and Kphy are pretty lengthy, 
we present them in the Appendix. These functions are plotted in Fig. [2] as these figures 
help us to get a better idea of how they behave as the parameters rjrei, qH* and r change. 
Besides their shapes, it is important to point out two issues. First, Umech is the same as the 
dispersion relation of the upper surface of a phase-separated binary low Reynolds number 
fluid with lower fluid viscosity rjs, thickness hs, and upper fluid viscosity r;^), thickness 
H* — hs- Second, in the limit rjrei = 1, Eq. fflTl) and flT^ become Eq. flTT]) and flT^ as we 
expected. 

Figure [2] shows Kmech and Kphy for rjrei 7 ^ 1- It can be seen that Kmech is always negative 
(stabilizing), and its magnitude increases with r. This is because the thickness of the lower 
(less viscous) region of the tissue increases and the damping is slowed down as r increases. 
On the other hand, Kphy is a non-monotonic function of r, as this term is originated from 
cell division/apoptosis, it cannot be simply explained by the viscous properties of the tissue. 
We will discuss more about the behavior of Kphy in the next section when the flow field in 
the tissue is examined. 

For fixed r and rjrei, Kphy has a minimum at qH* = 'E{r,r]rei)- Note that, unlike our first 
model, here S depends on f and r]reh Figure [3](a) shows S versus r for different choices of 
r]j.ei- It is clear that in our second model the minimum of ojphy also occurs at qH* ~ 0 ( 1 ). 
Figure [3](b) shows that when r and rj^ei are large, it is possible for Kphy to become positive. 
That is, it is possible for cell division/apoptosis to either slow down the relaxation toward 
the steady state, or even drive the tissue to an instability. This will be discused in details 
in Sec. IIII A1 

In this section rjrei and r are treated as independent variables. We will discuss the relation 
between rjrei and r in Sec. IIII B1 


III. RESULTS AND DISCUSSIONS 

A. Flow induced by cell apoptosis on the apical surface 

An immediate interesting question of our hydrodynamic theory of epithelium tissue is, 
besides fluidizing the tissue, how cell turnover inside the tissue affects tissue dynamics. This 
is considered in both our models. In the first model, we have assumed that the net cell 
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proliferation rate in a tissue is the only property that depends on the distance from the 
basal membrane, all other properties of the tissue are homogeneous. On the other hand, 
in the second model, regions with different net cell proliferation rate are assumed to have 
different viscosity. Interestingly, these two models predict very different tissue dynamics, 
revealing nontrivial effects of cell turnover on tissue dynamics. 

In our hrst model (Sec. IllAj) . the tissue is homogeneous except for the proliferative 
property. The linearized dynamics near the steady state shows that the steady state is 
always stable under small perturbations. The relaxation rate towards the homeostasis state 
contains two terms (Eq. (fTIJ]) h The hrst term oJmech (Eq. flTT]) l is driven by the surface tension 
of the apical surface, it is the same as the dispersion relation for the liquid-air surface of a 
low Reynolds number simple liquid. It always stabilizes the tissue and this rate increases 
with the wavenumber q (Fig. [T] (a)). 

The second term ijjphy (Eq. flT^ ) comes from cell turnover in the tissue. To gain further 
insight on the effect of this term to the dynamics of the tissue, in Fig. |4](a)(b) the held lines 
of Svy), the deviation of velocity held from the steady state in a tissue is plotted for a 
tissue with zero apical surface tension a for two diherent values of qH*. Since a = 0, the 
relaxation of the tissue is completely driven by the mechanism represented by ojphy It can 
be clearly seen that there is a how induced by cell apoptosis at the apical surface. As shown 
schematically in Fig. |4l(c), cell apoptosis close to the apical surface gives an ehective surface 
how. This is the mechanism that drives the how held shown in Fig. |l](a) and (b). Figure|ll(b) 
shows that when qH* is far from S (S ps 1.2), the how is restricted to region close to the 
apical surface. On the other hand. Fig. |4](a) shows that when qH* is closer to S, the how 
can penetrate deeper into the tissue. As q further decreases to qH* <C 1 the apical surface 
becomes nearly hat, and the induced horizontal how (not shown) is again small. Therefore 
the ehect of induced how is signihcant when qH* ^ S. This is rehected in Sec. Ill Al as the 
minimum of the magnitude of ojphy at qH* of order unity. The fact that ojphy has a minimum 
at qH* = S also indicates that although intuitively cell turnover helps to stabilize the tissue, 
how induced by apoptosis at apical surface slows down the relaxation of the tissue. Indeed, 
as can be seen from Fig. |4] (c), the induced how tends to push cells from regions where the 
tissue is thin to regions where the tissue is thick. Anyway, combining these two ehects, in 
this model ojphy is still negative for all q. That is, our hrst model predicts that cell turnover 
helps to stabilize the homeostasis state, although the stabilizing ehect is relatively weak 
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when qH* ~ 1. It is interesting to note that Fig. HJ^a) and Fig. HJ^b) also shows that vortices 
appears in the flow fleld when qH* > S, and they are absent when qH* < S. 

Fig. [5] shows the fleld lines of {Svx,Svz) predicted by our second model for a perturbed 
tissue with zero apical surface tension. Again the fleld lines for qH* > S look very different 
from the fleld lines for qH* < S. As qH* increases from qH* < S to qH* > S, vortices appear 
in the induced flow, this is similar to our first model (see Fig. 111(b)). Let he be the height of 
the center of the vortices. When qH* > S, Fig. [5] shows that for he > hs, Kphy decreases as 
f increases; however, when he < hs, Kphy increases as f increases. When he = hs, Kphy has 
a minimum (not shown in the figure). On the other hand, when qH* < S, Kphy increases 
monotonically with r. This shows that the shape of Kphy in Fig. [2](c) is related to both the 
flow induced by cell division/apoptosis and the thickness of the proliferative region. 

A very important result of our second model is that the sign of Kphy can be positive, i.e., 
it is possible that the overall effect of cell division/apoptosis is to hinder the relaxation of 
the tissue towards the steady state. This occurs when the following two conditions are both 
satisfied, (i) The viscosity rjs in the region h < hs is sufficiently small compared to rjo (the 
viscosity in the region > hs). (ii) The relative thickness of the proliferative region hs/H*, 
is sufficiently large. In Sec. IIII Bl we will discuss the relation between the distribution of 
proliferative cells and viscosity, there we will show that indeed in general rjrei = rjD/ris > 1. 
However, a tissue that satisfies condition (ii) is likely not a normal tissue, as a healthy 
mature epithelium tissue usually has a relatively thin layer of proliferative cells, a tissue 
with large hs/H* is more likely to be found during development or regeneration. 

Stability diagrams of our second model are shown in Fig. [6] (a) and (b). It is convenient 
to introduce a dimensionless parameter aq/{2ri£,r£)) to describe the relative importance of 
the apical surface tension and cell turnover in tissue dynamics. As can be seen by comparing 
Fig. |6] (a) and (b), tissues with small aq/ are relatively easy to become unstable, 
this is due to the relatively weak stabilizing effect from the apical surface tension. 

Intuitively cell turnover does not necessarily drive a tissue unstable. Indeed we know from 
our previous discussion that sometimes cell turnover helps to stabilize the homeostasis state, 
while sometimes cell turnover drives a steady state unstable. Therefore it is interesting to see 
when cell turnover stabilizes or destabilizes a steady state. Fig. |6] (c)-(f) shows how growth 
rate of small perturbations close to the steady state changes with our model parameters. 
It suggests that increasing rs and rjs helps to stabilize the steady state, as the peak of the 
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growth rate decreases; on the other hand, increasing ro and ? 7 d helps to destabilize the 
steady state, as the maximum growth rate of small perturbations increases. This further 
conhrms the intuition suggested by Fig. El^a) and El^b). 

Now we can summarize the role of cell turnover in the dynamics of a tissue close to a 
steady state. When the tissue height is perturbed from its steady state value, cell turnover 
drives the tissue toward the steady state. On the hand, cell apoptosis on the apical surface 
induces a flow in the tissue that brings more cells to thicker regions and pushes cells away 
from regions with smaller thickness. For a normal tissue, proliferative cells occupy a small 
region close to the basal surface, the overall effect of these mechanisms to tissue dynamics 
always helps to stabilize the tissue. On the other hand, for a tissue with a large proportion 
of proliferative cells, when the viscosity of the proliferative region of the tissue is sufhciently 
small compare to the viscosity of the non-proliferative region, the induced flow can destabilize 
the steady state, and the hrst instability occurs at wavelength 1 /g ~ H*. 


B. Cell division, cell apoptosis, and tissne viscosity 


Sec. IIIBI stated that the viscosity in a tissue depends on the local tissue composition. 
A derivation of the viscosity from a general hydrodynamic model of stratihed epithelium 
tissue is presented in the Appendix, where it is assumed that the stress relaxation process 
has negligible contribution from the turnover of junction proteins, and cell division/apoptosis 
produce force distribution that dominates the active stress inside a tissue. This approach 


different types of cells in a tissue 
can be expressed as 

T] = 


3 , 


generalizes the earlier theory of Ranft et ah by explicitly including the contribution from 


. Our analysis shows that the viscosity 77 in a tissue 


/i(To 


P 


{rsds - rDdD)^s + rodn 


( 16 ) 


where A 5 = ps/p, ps is the number density of proliferative cells and p is the total cell 
number density, p is the elastic shear modulus of the tissue, is the cell division rate and 
vd is the cell apoptosis rate, ds, do are characteristic strengths of the force dipoles during 
cell division and apoptosis, respectively, do characterizes the response of cell polarity to the 
stress in a tissue, cells in a tissue with small do are more likely to be polarized by mechanical 
stress. Equation flT^ suggests that viscosity in a tissue depends on the local proportion of 
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proliferative cells. This supports the argument we made in Sec. IIIBI 

Since in our model all the proliferative cells are localized in the basal {z < hs) region, 
A 5 = 1 for z < hs and A 5 = 0 for ^ > hs- The viscosity in the tissue can be expressed as 


7](z) = < 


J^sds ^ 




(17) 


From Eq. flTTl) . it is clear that the relative viscosity rjrei = Vd/vs and the relative proliferative 
rate r = ro/rs are not independent. In fact one can express rjrei as 

Vo 


'^rel — 


(18) 


where fjQ = ds/d-D > 0 is a. dimensionless parameter. 

Equation (fT7|) indicates that, for given fjQ, for a tissue to have large rjrei, the magnitude of 
f needs to be small. On the other hand, Eq. flT^ suggests that as r decreases, the proportion 
of proliferative region in the tissue also decreases. Putting these two results together. Fig. [7] 
shows rjrei versus r for a few choices of fjQ. For the contribution from Up^y to drive the tissue 
to an instability, the magnitude of fjo has to be large. Thus a tissue with small do is more 
likely to have a non-flat homeostasis state. 

Note that the viscosity shown in Eq. flTB]) is derived for a tissue whose cell adhesion protein 
turnover has negligible contribution to its dynamics. For a tissue in which the proliferative 
cells have small number of intercellular adhesion sites compared to terminal differentiated 
cells, rjrei could be greater than the magnitude predicted by Eq. ffTSj) . Therefore when both 
turnover of cell adhesion proteins and active force dipoles due to cell division/cell apoptosis 
contribute significantly to tissue dynamics, one still should treat r and rjrei as independent 
parameters. 

IV. CONCLUSION 


Although biological tissue exhibit high plasticity that allows remodeling, the homeostasis 
state is regulated such that its architecture is robust against intrinsic and external fluc- 
In this article we have studied how local tissue composition affects its 


tuations 12 


viscous properties, and how viscous properties in turn affects the dynamics of a tissue close 
to a steady state of a stratified epithelium. For a tissue with proliferative cells located 
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close to basal membrane, and differentiated cells located close to apical surface, the lower 
region of the tissue has a different viscosity compared to the upper region. Since cell di¬ 
vision/apoptosis inevitably induces a flow in this tissue, the position-dependent viscosity 
in the tissue has a profound effect on its dynamics. Our analysis shows that usually cell 
turnover is regulated to restore the steady state, but the flow induced by call apoptosis 
close to the apical surface induces a flow that hinder this relaxation process. As a result, 
perturbations with wavelengths comparable to the thickness of the tissue has the slowest 
relaxation rate toward the steady state. When the tissue has a thick proliferative region, 
and the viscosity of the proliferative region is signihcantly smaller than the other region, the 
steady state can even become unstable. 

Our model is quite general, thus our prediction should hold in general for stratihed 



result provides a possible route to trigger epithelial apical surface undulations. Although 
an instability due to the induced flow described by our model is unlikely to happen in 
a normal stratihed tissue, such mechanism could be important during embryogenesis and 
wound healing. Furthermore, when two different tissues are in contact 12|, the viscosity 
difference between these two tissues could also lead to interesting hydrodynamic instabilities. 
These possibilities will be explored in our future works. It is also worth mentioning that 
although in principle apical surface patterns can be a result of inhomogeneous distribution 
of growth factors in the tissue 3^, our model reveal that pure mechanical force is sufficient 


to induce apical surface patterns. 

To test our theory by in vitro experiments, one could try to increase the death rate 
of differentiated cells to increase r, and decrease the tissue apical surface tension to make 
the contribution from ujphy relatively more signihcant. In the case of mammalian olfactory 


epithelium, a higher death rate can be induced by unilateral olfactory bulbectomy [3^. The 
surface tension origins from the cell-cell adhesion [l, 2|, and it could be decreased by, for 
example the protease digestion procedure [40 1. 

Several modihcations can be made to make our model more general. For example, the 
solid stroma in our model can be easily replaced by a soft stroma. Another important 
improvement that we will make in the future work is to develop a model in whch cell 
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distribution is not pre-defined. It is important to check how homeostasis state and pattern 
formation is achieved in a stratihed epithelium with a more general hydrodynamic model. 
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work. 
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42| will have to be taken into account in this future 
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Appendix A: KpUy and Kmech 


1 . Kpfiy and Krnech 


The analytical form of Kphy and Kmech in Section 2.2 are 

{r]rei + lfi2qH* - sinhqH*) + {r]rei - Jis) 


K 


Kphy 


{r]rei + 1)^(1 + 2{qH*y + cosh2giP*) + {rjrei - hs) 

{r]rel + lf[2{qH*Y] + {r]rel - iffisiqH*, hs) 


- 1 , 


(Al) 


{r]rei + 1)41 + 2{qH*)‘^ + cosh2gTr*] + {r]rei - l) 2 p 2 (<?hf*, hs) 

where 

/ii =2qH*{l - 2hs[l + 2{qH*)%hs - 1 )^ 5 } + MH*Chs - 1) cosh{2qH*hs) - smh[2qH*{l - 2hs)] 
+ 2[1 + 2{qH*)^hl] smh[2qH*{l - hs)], 

(A2) 


P 2 =1 + 2{qH*)^l + 4(hs - l)h5[l + {qH*)\hs - l)hs]} + cosh[2gif*(l - 2hs)] 

- 2[1 + 2{qH*)^hl] cosh[2qH*{hs - 1)] - 2[1 + 2{qH*)\hs - 1)4 cosh{2qH*hs), 

(A3) 

P3 = 2{qH*)^{l + 2hs[-2 + (1 + 2{qH*)\hs - 1)4^5]} + 2(hs - 1)' cos\i{2qH*hs). (A4) 

Here hs = hs/H*. Although these expressions are complicated, they can be obtained from 
straightforward algebra. 
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Appendix B: viscosity in a tissue 


In this Appendix we first construct a general model for continuous self-renewing tissues 
that includes both spatial and cell lineage dynamics. By assuming that cell division and 
apoptosis dominates the stress relaxation in a tissue, we show that a tissue behaves as 
a viscoelastic material, and the viscosity in a tissue depends on local tissue composition. 
The parameters in the theory can be estimated from existing experimental data of model 
systems. The analysis in this Appendix follows closely the work of Ranft et ah ^ . The main 
difference is that in our model the cell lineage, which has been identihed as the fundamental 
unit of tissue and organ development [3^, is also taken into account. 

There are proliferative cells and terminal differentiated (TD) cells in a tissue. Typically, 
proliferative cells include stem cells and transit-amplifying (TA) cells differentiated from 
stem cells. For simplicity we ignore this internal conversion of proliferative cells, only mark 
that proliferative cells can undergo cell division, and the daughter cells have certain chance 
to become terminal differentiated cells. On the other hand, TD cells do not divide, they 


M, 


36 


41^ . we 


only undergo programmed cell death (apoptosis). Similar to other models 
neglect apoptosis of proliferative cells. 

Because of the coupling between the cortical network and the adhesion proteins, cells in 
a tissue are tightly connected { 2 , 5|, and diffusion of cells inside the tissue can be neglected 
44l |. Denote ps as the number density of proliferative cells, po as the number density of TD 
cells. Taking into account cell division/apoptosis, and advection of cells by the flow inside 
the tissue, the continuity equations for ps and po can be written as 


dtps + di{vips) = rs{2ps - l)ps, 


(Bl) 


dtpD + diiyipo) = 2rs(l - ps)ps - tdPd, (B2) 

where Einstein summation convention is used to simplify the notation, rs and td are the 
division rate of proliferative cells and apoptosis rate of TD cells, respectively, ps is the 
self-renewal probability of proliferative cells, and Vi {i = x, y, z) is the Ath component of the 
flow field. 

The total cell density at position r is p(r, t) = ps{r, t) -|- poi'^, t)- The fraction of the pro¬ 
liferative cells at r is As{T,t) = ps{r,t)/p{r,t). From Eq. flBlDflB2p . the evolution equations 
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for p and A 5 are 

Dtp = p[{rs+ rj:))Ks - re - divi], (B3) 

Dtks = [^rsPs -rs + ro- {rg + rD)^s\^s, (B4) 

where Dt = dt + vidi is the material derivative. 

We assume that, in the absence of active processes such as cell division and cell apoptosis, 
a tissue behaves as a linear elastic material, thus there is a linear relation between the change 
of elastic stress Acr^® and the change of elastic strain Auik, Aafp. = {K — 2u/S)AunSik + 
2pAuik, where K is the compressional modulus, and p is the shear modulus I45II . Besides 
elastic stress, cell division and cell apoptosis also exert forces to the tissue Because 

these processes exert no net force and torque on the tissue l48|, we model them as 
symmetric force dipoles [?, ^|. For a cell division/apoptosis event occurring at r = tq, 

the change of stress is related to the active force dipole dik through = —dikS{r — Tq), 


where 5(r — tq) is the Dirac delta function, and superscript ‘D 
active cell division/apoptosis events. Since inertia is negligible 
gives 

0 , 

where atk = afk + is the total stress in the tissue. 


denotes stress created by 


49[ |. force balance condition 

(B5) 


Due to cell division and apoptosis, a unique reference state for the strain does not exist. 


3 - 


However, the difference of strain between subsequent states still can be dehned [711. Let the 
frame of reference flow and rotate with local flow and vortex, the evolution equation for the 
stress can be written as 


D/cTjfc = (A - ^p)vii6ik + 2pvik + D^afk, 


(B 6 ) 


where Dtatk = Dtatk + Duaik + Dkiau is the Jaumann derivative. Here Vik = {diVk + dkVi)/2 
is the strain rate tensor, and = {diVk — dkVi)/2 is the vorticity tensor. 

It is convenient to decompose the total stress tensor into isotropic and traceless parts, i.e., 
o'ik = —pdik + dik, where p = —au/S is the tissue pressure, and dik is the traceless (deviator) 
part of total stress tensor. Physically the isotropic part of the stress tensor is related to the 
change of local tissue volume, and dik is related to the distortion of local tissue shape [45 1. 
The same decomposition can also be applied to active stress tensor and strain rate tensor, 
i.e. = —p^dik + and Vik = (l/3)n;Aifc + Dfc- Now the stress evolution Eq. flB 6 p can be 
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expressed as 


Dtp = -Kvu + Dtp-', 


(B7) 


and 




ik * 


(B8) 


To complete our model, we still need the constitutive equations for Dtp^ and D^a^. 
Active stress is induced by cell division and cell apoptosis, we first consider the force dipole 
created by cell division. Let unit vector /j denote the direction of the cell division mitotic axis 
(Fig. |8]), follow 1^, the force dipole produced by a cell division event can be approximated 
by 

(B9) 




^ik 

Here ds characterizes the strength of the dipole. Typically > 0 because cell division 
makes tissue grow. The magnitude of ds can be deduced from experiments. Since different 
dividing cells may have different mitotic axis, the macroscopic properties is related to the 
average of dik in a small region. 


{dfk) — ^dtsiklk) — ds{5ik + dik), 


(BIO) 


where qik = (Slj/fc — dik) is the nematic order tensor in liquid crystal literatures. In general, 


tissue growth and cell division/apoptosis processes are intrinsically anisotropic 
and the orientation of the cell motitic axis depends on local tissue distortion [7|, 
linear order, this relation can be expressed as 

_n ^ik 

Qik - Qik = -, 

<70 


i 
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56 1. To 


(Bll) 


where do > 0 is a constant, and is the intrinsic anisotropic tensor, which should satisfy 
q^i = 0. Denote the eigenvalues of qik — q^k as q^^\ it is convenient to write 


= 3(cos^ 6i) — 1, 


(B12) 


where 6i is the direction angle of mitotic axis relative to the Tth principle axis of qik — q^k- 
From Eq. (IBllll we have 




(k) 


(B13) 


k=l 
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where Si = (3(cos^6'i) — l)/2 is the three-dimensional order parameter relative to i-th 
principle axis of qik — and cr’^*^’s are the i-th. eigenvalue of aik- From this relation, one 
can estimate ao from experimental measurements 56| . 

Cell apoptosis, a key mechanism of tissue morphogenesis [^, , also exert force on the 

surrounding environment (Fig. [8]). During apoptosis, a cell rapidly develops an actomyosin 
ring around its periphery and signal to neighboring cells to induce “purse-string”-like con¬ 


tractility in the neighboring cells 


59| . This “purse-string” contraction depends on myosin 




activity, which can be elevated by applying mechanical forces to the tissues 60|. Therefore, 
in general the force dipole exerted by cell apoptosis is anisotropic and dependent on the 
local stress. 

{<r> = + do f?! + , (B14) 

the traceless part of this force dipole has an intrinsic contribution proportional to and 
a contribution induced by local stress, dn and are constants characterizing the strength 
of the isotropic and the traceless part of the dipole. Typically do < 0 since cell apoptosis 
often leads to the decrease of tissue volume. On the other hand, there is no simple intuitive 
argument to determine the the sign of d^. 

Substitute Eq. flB10j)flB14|) into Eq. flB7j)flB8D leads to 


DtP = -Kvii + p[{rsds - rDdD)J^s + 


(B15) 


(1 -h TD^)dik = 2r]Vik + dlf,, (B16) 

where = Tp[{rsds — rj:)dD)As + rDdD]Qik is the intrinsic active stress due to cell intrinsic 
polarity, p = rp, and the time scale r that characterizes the active stress induced by cell 
proliferation/apoptosis is given by 


r ^ =—[{rsds-rDdD)As + rDdD]- 
o'o 


(BIT) 


From Eq. flB16j) and Eq. flB17j) . we make the following remarks. 


For a tissue with no intrinsic cell polarity, = 0. The stress evolution Eq. flBlbj) 
reduces to the Maxwell model of a visco-elastic material with stress relaxation time r 
and viscosity p. At time scale greater than r, the tissue behaves like a fluid j^. 


difc = 2pvik. 


(B18) 
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It is important to note that the viscosity of the tissue depends on p and A 5 . This 
supports the intuitive argument that we presented in the main text that the tissue 
viscosity depends on the local density of proliferative cells and TD cells. 


From Eq. flBlTp a sufficient condition for the tissue viscosity to be positive is do > 0 


This is consistent with what has been observed in the experiments 8l4l0|. 


In general is not zero. In this case, in the long time limit the stress of the tissue 
satisfies 

dik = 2r]Vik + d/fc. (B19) 

The additional term on the right-hand side of the above equation suggests that a 
tissue in a flow-free state is in general under stress. This is the active stress from cells 
with intrinsic polarity. It is important for tissue morphogenesis, for example, in the 


imaginal disk of Drosophila flisai. 

At time scale large compare to r, the total stress tensor in a tissue reduces to 

(Tik = -pdik + 27]Vik + d/fc. 


(B20) 


Hence a tissue in this limit behaves as a viscous fluid coupled to a scalar field A 5 and an order 
parameter which measures the intrinsic anisotropy of cell orientation. Moreover, Eq. flB20ll . 
Eq. (lB3tlB5[l and Eq. (IB 151) give a close set of equations describing tissue and cell lineage 
dynamics. Note that the constitutive relation we used in the main text describes a tissue 
with negligible intrinsic active stress. 

By measuring traction stress exerted W the dividing Dictostelium cell on the flexible 
poly-acrylamide gel, Tanimoto and Sano 5^ found that the typical strength of force dipole 
during cell division is about 1.2 x 10“^^ N-m and the orientation of the force dipole coincides 
with cell division axis 5^. This suggests that ds in Eq. flB9ll can be measured. Although 
Dictostelium does not form tissues, one can use this experimental measurement to get a 
rough idea of the order of magnitude of ds as ds ~ 4 X 10 N ■ m. 

Experiments similar to 5^ can be applied to measure (Tq. In this experiment, Marel et al. 
measured the average cell division orientation as a function of local strain rate in monolayer 
Madin-Darby canine kidney (MDCK) cells and concluded that 


2(cosV) - 1 = |(Ai - A 2 ), 


(B21) 
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where Ai and A 2 (Ai > A 2 ) are the eigenvalues of strain rate tensor and 0 is the angle 
between cell division axis and the principle axis of the eigenvalue Ai. By assuming the 
monolayer tissue behaves as a Newtonian fluid, we can use the two-dimensional form of 
Eq. dmsD , 


(B22) 


2(cos2 0)-1^^(Ai-A2), 
<^0 


where fj is the 2-dimensional viscosity. Using typical tissue viscosity rj = 10® Pa • s [10| and 
MDCK monolayer height h = 2 fim a comparison between Eq. (IB21j) and Eq. flB22j) 
gives do = 5 Pa • m. Thus do = ao/h = 2.5 x 10® Pa. 

For amnioserosa tissue, the magnitude of do can be inferred from the apoptosis rate 
measured from dorsal closure process of Drosopila. The measurement made by Toyama et 


-5 ^-1 


ah gives = 5.5 x 10 ® s 
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,Q- 


Since during dorsal closure, amnioserosa cell only 


undergo cell apoptosis, Eq. flB17D is reduced to r ^ = prndu/aQ. Use tissue heigh 


h = 5/im 


63j |. cell density p = 2.2 x 10 ® pm ^ [57|, and the relaxation time r = 100 s [64|. If we 
use the magnitude of (Tq measured from MDCK cells, we find do — 2.1 x 10“^ N ■ m for 


amnioserosa tissue. In principle, combining the measurement by 
like 5^ for amnioserosa tissue give us a better estimate for its 


57l | with an experiment 


In summary, in this Appendix we have presented a theoretical model to relate the viscosity 
of a tissue to a few measurable parameters. The key parameters d^, d/), and do of a tissue can 
be measured from different experiments. Therefore it is possible to infer the hydrodynamic 
properties of a tissue from existing experimental techniques. 
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Atotai is the number of amnioserosa cells just before dorsal closure. Typically Atotai ~ 100. 
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Scaled relaxation rate 



qH* 


FIG. 1. 2r]UJmech/(^Q (blue curve) and ujphy/kp{H*) (red curve). 0Jmech{<l) decreases monotonically 
towards (Tq/2r] as qH* ^ 1, \ujphy\ has a minimum at qH* ~ 1. 
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FIG. 2. (a) Kmech as a function of f and qH* for ij^ei = 5 . (b) K^ech versus r for qH* = 2, 

Vrei = 5 (solid blue curve), qH* = 0.6, = 5 (solid red curve), qH* = 2, rjrei = 4 (dashed 

blue curve), and qH* = 0.6, tj^ei = 4 (dashed red curve), (c) Kp^y as a function of f and qH* for 
rjrel = 5 . (d) Kphy/ro versus r for qH* = 2, rj^ei = 5 (solid blue curve), qH* = 0.6, rj^ei = 5 (solid 
red curve), qH* = 2, = 4 (dashed blue curve), and qH* = 0.6, rj^ei = 4 (dashed red curve). 
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FIG. 3. (a) H as a function of r for r]rei = 0.4 (blue), 1 
the least negative value of Kphy for given f, rjrei, i-e., 
shows where Kphy{qH* = H,f, r/^ez) is zero. 


(red), and 5 (orange), (b) This plot shows 
KphyiqH* = E,f,r]rei)- The dotted curve 

























FIG. 4. (a) and (b): Field lines of {6vx, 6vy) calculated from our first model for a perturbed tissue 

with tissue apical surface tension u = 0. It can be seen that the topology of the flow field depends 
on the magnitude of qH*. qH* = 0.8 in (a), qH* = 2 in (b). The red arrows are the stream lines, 
the dashed black lines are the apical surface of the perturbed tissue, and the background gray-level 
indicates the magnitude of 6v. To generate the flow field, we have chosen kp{H*) = —1, H* = 2, 
and 6H = O.Olcos(gx). (c) Cell death at apical surface induces a flow held that is indicated by 
the empty arrows. Close to the apical surface this how held has a horizontal component indicated 
by the blue arrows. The solid curve is the apical surface at the present moment, and the dashed 
curve indicates the apical surface at a later time. 
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FIG. 5. Field lines of {6vx,6vy) in the second model for cr = 0, and Kphy versus f for qH* = 0.8 
(< H for all f) and qH* = 2.0 (> S for all f). For qH* = 0.8, Kphy increases monotonically with 
r. For qH* = 2, there are vortices in the induced flow. Kphy decreases with f when he > hs, Kphy 
increases with f when he < hs- hs is indicated by the green lines, and he is the z-position of the 
vortex center. To generate the flow fields, we have chosen rjeei = 5, kp{H*) = —1, H* = 2, and 
6H = 0.01 cos {qx). 30 
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FIG. 6. (a)(b) Gray scale indicates th relaxation rate of a perturbed tissue scaled hy ru for 

qH* = H. (a) aq/{2riDrD) = 0.4, and (b) aq/{2rioi'D) = 0.05. (c)-(f) relaxation rate (in units 
of 1/day) of a tissue as a function of the wavelength of perturbation. Solid curves; r 5 = 1 d“^, 
ri) = 30 d“^, ris = 10 MPa s, rjo = 50 MPa s, cr = 1 m Nm“^, and hs = 30 //m (the corresponding 
magnitude of H* is 31 //m). The dashed lines give the relaxation rate versus A for the same 
parameters except rs = 2 d“^ in (c), ro = 43 d“^ in (d), r]s = 20MPa s in (e), and ri£) = SOMPa s 
in (f). 
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FIG. 7. r versus r]rei for fjQ = 30 (blue curve), 90 (red curve), and 150 (yellow curve). The short- 
dashed curve is the stability boundary for a tissue with aq/{2r]Dro) = 0.2 and qH* = H. The 
stability boundary is determined from the zero of the relaxation rate for perturbations with the 
most dangerous wavelength. 
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FIG. 8. Cell division and apoptosis as a source of localized force inside tissue, (a) Cell division 
can produce equal but opposite forces along the mitotic axis (dashed arrow), (b) Cell apoptosis 
induce a purse-string (dashed line) which shrinks the cell surface and prevents gap formation. Both 
processes create a force dipole acting on the surrounding environment (blue cells). 


33 

















